Fisheries genomics of snapper (Chrysophrys auratus) along the west Australian coast

Abstract The efficacy of fisheries management strategies depends on stock assessment and management actions being carried out at appropriate spatial scales. This requires understanding of spatial and temporal population structure and connectivity, which is challenging in weakly structured and highly connected marine populations. We carried out a population genomics study of the heavily exploited snapper (Chrysophrys auratus) along ~2600 km of the Australian coastline, with a focus on Western Australia (WA). We used 10,903 filtered SNPs in 341 individuals from eight sampling locations to characterize population structure and connectivity in snapper across WA and to assess if current spatial scales of stock assessment and management agree with evidence from population genomics. Our dataset also enabled us to investigate temporal stability in population structure as well as connectivity between WA and its nearest, eastern jurisdictional neighbour. As expected for a species influenced by the extensive ocean boundary current in the region, low genetic differentiation and high connectivity were uncovered across WA. However, we did detect strong isolation by distance and genetic discontinuities in the mid‐west and south‐east. The discontinuities correlate with boundaries between biogeographic regions, influenced by on‐shelf oceanography, and the sites of important spawning aggregations. We also detected temporal instability in genetic structure at one of our sites, possibly due to interannual variability in recruitment in adjacent regions. Our results partly contrast with the current spatial management of snapper in WA, indicating the likely benefits of a review. This study supports the value of population genomic surveys in informing the management of weakly structured and wide‐ranging marine fishery resources.


| INTRODUC TI ON
Marine ecosystems are one of the last on the planet in which wild populations are heavily exploited for human consumption. Not only do aquatic animals constitute ~17% of the world's meat consumption, but their capture supports the livelihoods of an estimated 10-12% of the global population and many species hold cultural significance (FAO, 2020). Although more than half of all assessed fish stocks were likely depleted by the mid-1970s, the proportions of depleted stocks recently dropped to below one quarter in regions practising intensive evidence-based fisheries management (FAO, 2020;Hilborn et al., 2020). Indeed, the intensity of management efforts and knowledge of species biology strongly correlate with fisheries sustainability and therefore realized fisheries yield (Hilborn et al., 2020;Nilsson et al., 2019). Knowledge of biological stock structure (i.e. population structure), defined as the number and spatial extent of distinct populations within a species' range, as well as the levels of connectivity between such stocks, is essential for the performance of fishery assessments and sustainable management (Cadrin, 2020). This is because distinct stocks may display different dynamics (e.g. recruitment, growth rates), demographics (e.g. sex ratios, abundance) and genetics (e.g. diversity, environmental adaptations), causing unique responses to fishing and environmental pressures . The consequences of not incorporating accurate information on stock structure and connectivity in assessment and management have been widely documented (Fu & Fanning, 2004;Kell et al., 2009;Kerr et al., 2014;Smedbol & Stephenson, 2001;Sterner, 2007;Ying et al., 2011). For example, consideration of the complex stock structure of Bristol Bay sockeye salmon (Oncorhynchus nekra) has been important for maintaining fishery productivity, while depletions have occurred in other places where spatial structure has not been managed effectively (Hilborn et al., 2003;Schindler et al., 2010).
Key biological traits that shape population structure and connectivity in marine species include: habitat specialization, determined by breeding site preferences, dietary requirements and abiotic tolerance limits (Cowen, 2002;Cowen & Sponaugle, 2009); and dispersal potential (Bohonak, 1999), determined by reproductive strategy (e.g. broadcast spawning and brooding), the duration of any pelagic stages and the movement behaviour of both larval and postlarval stages. Due to the fluidity of the marine environment and processes like currents, waves and tides, connectivity and population homogeneity can occur over large spatial scales, particularly in species with long pelagic larval stages (PLS; Cowen & Sponaugle, 2009;Selkoe et al., 2008). For example, in the broadcast spawning and highly migratory yellowfin tuna (Thunnus albacares), population homogeneity likely occurs at spatial scales as extensive as entire ocean basins Pecoraro et al., 2018). Nevertheless, marine environments are dynamic and heterogeneous and therefore contain elements that can disrupt connectivity and act as barriers to dispersal, including meanders and eddies, fronts, irregular coastline topology, habitat heterogeneity and countercurrents (Cowen, 2002;Cowen & Sponaugle, 2009). For example, Taillebois et al. (2017) detected significant population structuring in the broadcast spawning marine teleost, the black-spotted croaker (Protonibea diacanthus), across topographically complex coastline in northern Australia at spatial scales of only 100 s of kilometres.
Methods for assessing stock structure and connectivity in exploted marine species include otolith microchemistry, isotope analysis, mark-recapture and DNA analyses. Advances in DNA sequencing technologies have increased the value of DNA-based methods in stock delineation studies, and it is now economical to generate large datasets of 1000 s of DNA markers (i.e. those used in genomic studies), vastly contrasting with previously utilized datasets of 10 s of markers (i.e. those used in microsatellite studies). Large datasets of DNA markers have greater resolving power to uncover F I G U R E 1 Map of the study region showing the eight sampling sites, including the comparative site in South Australia (SA), the three current management areas for snapper in Western Australia (WA): Shark Bay Oceanic (red), West Coast (blue) and South Coast (green), as well as the westernmost management area in SA: Spencer Gulf/West Coast (purple). Management area boundaries are indicated with dotted black lines. The inset plots summarize the lengths (in mm) and ages (in years) of sampled snapper (excluding the temporal samples from CS and ESP). Sampling sites in WA: GAS, Gascoyne; KAL, Kalbarri; LAN, Lancelin; CS, Cockburn Sound; BUS, Busselton; ALB, Albany; ESP, Esperance. Sampling site in SA: CED, Ceduna the often subtle patterns of population structure and connectivity characteristic of marine species (Bernatchez et al., 2017;Grummer et al., 2019), occurring due to their high abundances and dispersal abilities.
Our focus here is on snapper (Chrysophrys auratus), a large, long-lived demersal sparid distributed in subtropical and temperate coastal waters of Australia and New Zealand (Gomon et al., 2008). In Western Australia (WA), snapper supports highly valuable commercial and recreational fisheries from the Shark Bay region to the border of South Australia (SA; Figure 1). In 2017/18, snapper production by the WA commercial sector was valued at ~1.7 million AUD (Steven et al., 2021), while snapper constituted ~4% of catch by the WA recreational sector , an industry that generates ~2.4 billion AUD per year (McLeod & Lindner, 2018). During 2017/2018, ~242 tonnes of snapper were landed in WA state waters across all fishing sectors (of which ~50% was commercial catch), with the largest catches landed off the west coast (Gaughan & Santoro, 2020;Ryan et al., 2019). Due to its fishery importance, its economic and social importance, and its inherent vulnerability to fishing as a result of its biology (e.g. long-lived, aggregate spawner), snapper is used as an indicator species for the inshore suite of demersal scalefish resources across its WA range (Newman et al., 2018). In WA, snapper is assessed and managed as three separate stocks (not including the three stocks in the inner gulfs of Shark Bay not addressed in this study)-the Shark Bay Oceanic, West Coast and South Coast stocks ( Figure 1). The Shark Bay Oceanic and West Coast stocks are currently considered depleted, with management in place for recovery, while the South Coast stock is considered sustainable (Fairclough et al., 2021;Newman et al., 2022). These stocks and their spatial boundaries are based on bioregions defined primarily from information on environmental characteristics rather than from knowledge of species-specific population structure and stock connectivity (Newman et al., 2022).
The particular reproductive characteristics and movement behaviours of snapper as well as the dominant patterns of oceanographic circulation along the WA coast are likely important factors for shaping stock structure and connectivity in the species. The dispersal potential of snapper is expected to be high, particularly during the pelagic and subadult stages (Fairclough et al., 2013;Parsons et al., 2014). Snapper are multiple batch, broadcast spawners, with eggs hatching after 1-2 days and the pelagic larval stage lasting for ~17-30 days (Fowler & Jennings, 2003;Francis, 1994). When eggs and larvae occur in open waters, they are expected to be influenced by the Leeuwin Current (LC) and the seasonal, wind-driven Capes and Cresswell Currents (see Figure 1, Akhir et al., 2020;Cresswell & Golding, 1980;Cresswell & Peterson, 1993;Pearce & Pattiaratchi, 1999). After settlement, juvenile snapper migrate from spawning areas to sheltered inshore habitats where they reside until ~2 years of age. During the subadult stage, snapper are known to migrate distances of more than 1000 km to the coastal regions where they become resident (Fowler et al., 2005;Hamer et al., 2011;Wakefield et al., 2011).
Other life-history and oceanographic characteristics could act to limit connectivity and dispersal in snapper in WA, leading to population structure. Although snapper spawn along the entire WA coastline, a number of embayments where large spawning aggregations occur are particularly important (Nahas et al., 2003;Wakefield, 2010;Wakefield et al., 2015). These include coastal embayments of Shark Bay, Perth (Cockburn Sound, Warnbro Sound, Owen Anchorage) and Albany (King George Sound; Figure 1). These environments are largely protected from oceanographic processes in open shelf waters and also often feature local circulation patterns that act to retain eggs and larvae (Steedman & Craig, 1983;Wakefield et al., 2011). With respect to older life-stages, although adult snapper are known to travel distances of hundreds of kilometres or more, most show more limited movements of spatial scales of <100 km (Crisafulli et al., 2019;Moran et al., 2003;Sumpton et al., 2003).
Current knowledge of population structure and connectivity in snapper in WA is based on tagging, mark-recapture, otolith microchemistry and microsatellite analyses. Recent mark-recapture work on the lower west coast of WA indicated that adult movement is largely limited to spatial scales of <20 km, although this estimate is probably influenced by spatial variability in fishing effort (Crisafulli et al., 2019). Microsatellite DNA analyses have suggested that snapper in WA are characterized by an isolation by distance pattern of population structure, rather than genetically distinct subpopulations (Gardner & Chaplin, 2011). Otolith microchemistry studies have reported different patterns of mixing during different life stages as well as a range of patterns of spatial differentiation using different chemical tags (Edmonds et al., 1999;Fairclough et al., 2013). Additionally, the latter study (i.e. Fairclough et al., 2013) demonstrated that adults in any one location are derived from multiple nearshore nursery environments.
Therefore, a degree of uncertainty still remains around where appropriate stock boundaries should be drawn for management and assessment purposes and whether current boundaries are appropriate. In addition, declines in snapper catches and stock depletions in parts of WA point to the need for more information about population structure and connectivity across the state (Fowler et al., 2021).
In this study, we use genome-wide polymorphism data and population genomic analyses to characterize population genetic structure and connectivity in snapper across its WA range. Our dataset also enables us to assess temporal stability in population structure as well as connectivity in snapper between WA and its nearest, eastern jurisdictional neighbour (SA). Given the species' life-history traits and the oceanographic setting of WA's coast, we predict population differentiation in snapper to be influenced by broad patterns of on-shelf oceanographic circulation. Our secondary goal is to determine whether current spatial scales of assessment and management reflect the biological units (i.e. stocks) identified with genomics, as well as their spatial boundaries. Our third goal is to demonstrate to fisheries scientists and managers the value of genomic datasets in clarifying stock structure and connectivity in a species for which this was previously difficult (i.e. highly abundant and dispersive species using smaller genetic datasets).

| Sampling and associated biological information
Muscle or fin-clip samples were obtained from 290 snapper from eight locations between the Gascoyne, on the central west coast of Western Australia (WA), and Ceduna, on the west coast of South Australia (SA; see Figure 1, Table 1 and Table S1). Tissue samples were obtained from fish landed by recreational or commercial fish- Our localities also cover the management areas used for snapper in WA outside of the inner gulfs of Shark Bay (not addressed here) -Shark Bay Oceanic, West Coast and South Coast (see Figure 1).
Tissues were preserved in 100% ethanol and stored at −20°C until DNA extraction. Where possible, biological data including total and fork length, age, sex and reproductive stage were obtained for each sampled individual (see Table 1 and Table S1).
Additional tissue samples were obtained from snapper landed off Esperance (n = 28) and in Cockburn Sound (n = 30) in 2010 and 2014, respectively, to investigate the temporal stability of trends in population structure. Although the time periods covered preclude investigation across multiple generations, they are adequate for assessing interpopulation migration and population range shifts as a result of temporal variation in recruitment or population density (information important in fisheries management). The Cockburn Sound samples were collected by DPIRD WA staff as part of a fisheries independent survey, while the Esperance samples were collected by DPIRD WA staff from commercial catch for a former population genetics project (see Gardner et al., 2014Gardner et al., , 2017Gardner & Chaplin, 2011).

| DNA extraction, genomic library preparation and sequencing
Genomic DNA was isolated from each sample using a modified salting-out protocol (Sunnucks & Hales, 1996). Quality control of DNA extracts was carried out with NanoDrop and gel electrophoresis. Extracts were quantified with Qubit and diluted to ~10-

| Bioinformatics
Raw sequence reads were processed to generate a high-quality SNP dataset using similar bioinformatic procedures as detailed elsewhere (Sandoval-Castillo et al., 2018), but with the assistance of a reference genome for snapper. Specifically, the quality of raw sequence data was checked using FastQC before being demultiplexed with the process_radtags module from STACKS 2.0 (Catchen et al., 2013).
Barcodes, restriction sites and RAD tags were then trimmed from sequence reads using TRIMMOMATIC (Bolger et al., 2014). Trimmed sequence reads were then aligned to a high-quality snapper reference genome (Catanach et al., 2019) using BOWTIE 2 (Langmead & Salzberg, 2012). The SNPs were subsequently called using BCFTOOLS (Narasimhan et al., 2016). The resulting dataset was initially filtered using VCFTOOLS (Danecek et al., 2011) to retain only bi-allelic SNPs present in at least 80% of individuals in all populations with a minimum minor allele frequency of 0.03. Also using VCFTOOLS, further filtering was carried out to remove indels, individuals with more than 20% missing data, SNPs with low and extremely high coverage, SNPs with low mapping quality, SNPs not in Hardy-Weinberg Equilibrium (HWE) and physically linked SNPs.

| Categorizing putatively neutral SNPs
The Bayesian method in BAYESCAN 3.0 (Foll & Gaggiotti, 2008) was used to identify candidate SNPs putatively under selection. The software was run with 20 pilot runs, each with 5000 iterations, followed by 100,000 iterations with a burn-in length of 50,000 iterations. The outlier SNPs were identified using a 5% false discovery rate with a prior odd of 10 and were subsequently removed from the dataset to produce a putatively neutral one. The study of the role of natural selection on snapper populations using outlier SNPs and genomic regions associated with environmental variation is the topic of a separate and ongoing investigation (Brauer et al. unpublished). Missing genotypes (~0.6% of data matrix) were assigned with the most common genotype at that locus. A discriminant analysis of principal components (DAPC) was also conducted using the R package ADEGENET 2.1.5 (Jombart, 2008;Jombart et al., 2010)

| Isolation by coastal distance
We tested for a signal of isolation by distance (IBD) across the sampling range by assessing the relationship between coastal distance with a Mantel test (9999 permutations) in GENALEX 6.5 (Peakall & Smouse, 2012). Distances between sampling locations were estimated as the shortest distance between sites following the coastline using the viamaris function in MELFUR 0.9 (https://github.com/pygmy perch/ melfuR). Mean effective dispersal distance was then estimated from the slope of the IBD relationship using the theoretical model of Kinlan and Gaines (2003), an extension of Palumbi (2003) ies management (Palumbi, 2003). Of further relevance to fisheries management is that IBD slopes, unlike F ST , are not affected by rare dispersal events but strongly reflect dispersal over the most proximate generations and therefore over ecologically significant time scales (Rousset, 1997).

| Connectivity barriers
We tested for putative barriers to dispersal using the piecewise regression approach implemented in Robinet et al. (2020) to identify regions where gene flow deviates from an IBD pattern. Briefly, the individual ancestry proportions generated by ADMIXTURE for a K of 2 (presented in Figure 2c) were used to estimate a mean SA ancestry

| Local-scale structure
To investigate local-scale patterns of gene flow and further explore the impact of coastal distance on genetic structure, we assessed the genetic similarity between individuals at increasing geographic distances with spatial autocorrelation analyses in GENALEX 6.5 (Peakall & Smouse, 2012;Smouse & Peakall, 1999). Spatial autocorrelation analysis can provide higher resolution information on current patterns of gene flow than evolutionary estimators like F ST (Peakall et al., 2003). Specifically, the analysis can uncover IBD signals over smaller scales and demarcate ecologically important genetic patches. The extent of nonrandom and positive spatial autocorrelation can be inferred from the first x-intercept in the correlogram (also referred to as genetic patch size) if a significant correlation coefficient (r) occurs in at least one distance class (Smouse & Peakall, 1999;Sokal & Wartenberg, 1983).

| Genetic diversity, population differentiation and clustering analyses
Levels of genetic diversity were very similar across sites ( Table 1). The PCA indicated the presence of three geographically distinct groups across the study region, referred to herein as the midwest (GAS, KAL, LAN), south-west (CS, BUS, ALB) and south-coast (CED; Figure 2b). The ESP sample did not cluster predominately with any one group but spanned across the south-west and southcoast. Differentiation between the groups varied, with greater distinction observed between the south-coast and the south-west and mid-west groups than between the mid-west and south-west groups. Of the two WA groups, the mid-west was more tightly clustered than the south-west, suggesting greater homogeneity in the former. The DAPC results were consistent with these PCA outcomes ( Figure S1).

Expected heterozygosity (H
The ADMIXTURE analysis suggested that the most probable number of genetic clusters in the dataset was two (i.e. K = 2). These clusters loosely corresponded to the samples between GAS and ALB Comparing the two temporal CS and ESP samples, membership to the three groups identified with ADMIXTURE indicated temporal stability in genetic structure at the former location only (Figure 3).
Compared with the ESP sample, ESP10 had higher membership to the south-coast group (Wilcoxon test: p-value < 0.0005) and lower membership to the mid-west and south-west groups (Wilcoxon test: p-value = 0.001 and 0.02, respectively).

| Isolation by coastal distance
We detected highly significant IBD across the sampling range (r = 0.86, p-value <0.001; Figure 4). This analysis indicated that spatial distance explains 74% of the variation in linearized F ST, accounting for a substantial amount of the population genetic differentiation inferred across the sampling range. The mean effective dispersal distance per generation was estimated as 400 km from an IBD slope of 4E-06 (i.e. F ST = 0.004 per 1000 km). This suggests that demographic connectivity between locations separated by distances greater than 400 km may be limited. Our estimate of 400 km is within the range of those made for other marine fish (Kinlan & Gaines, 2003).

| Barriers to connectivity
The These results suggest that the genetic structure observed along  Figure 2c) across the WA samples (±standard error), determined from the ADMIXTURE results for K = 2, as a function of distance from ESP, and (b) results of the piecewise regression analysis showing residual standard error as a function of distance from ESP the WA coast is not merely due to IBD but is also a reflection of a connectivity barrier, an interpretation consistent with the results of the clustering analyses.

| Local-scale structure
The spatial autocorrelation analyses indicated that across the sampling range, positive genetic autocorrelation occurs between individuals sampled at the same site (Figure 6c), a result suggestive of recruitment to the local subpopulation. Within site spatial autocorrelation was highest at GAS and KAL in the mid-west and at ALB in the south-west (Figure 6c). Variation around samplespecific r values was highest for KAL, ALB and ESP. These samples were the result of the greatest number of fishing days that were separated by the greatest stretches of time (Table S1). This could perhaps indicate that cohesion among different groups of genetically alike individuals may occur in snapper across this region, but such a possibility would need to be verified with larger samples than those available here.
Our spatial autocorrelation correlograms indicated that localscale genetic structure occurs within the south-west but not the mid-west (Figure 6a

| DISCUSS ION
This is the first study to investigate population structure and connectivity in the fishery important Australian snapper using a population genomics approach. We detected spatially variable population genetic structure and levels of connectivity in snapper across ~3000 km of coastline. Although our results indicated that snapper across Western Australia (WA) are characterized by weak genetic structure, we found evidence for stock discontinuities, dispersal limits and local-scale structure. Our results improve upon prior understanding of population structure and connectivity in snapper across its WA range based on tagging (Crisafulli et al., 2019;Wakefield et al., 2011), otolith microchemistry (Edmonds et al., 1999;Fairclough et al., 2013) and microsatellite analyses (Gardner & Chaplin, 2011). They also partly contrast with the spatial scales of assessment and management currently used for snapper in WA (see Fairclough et al., 2021), indicating that a review may be beneficial.
Our study demonstrates the value of population genomic surveys to inform the management of weakly structured and wide-ranging marine fishery resources.

F I G U R E 6
Spatial autocorrelation analyses of snapper samples from the (a) mid-west, (b) south-west and (c) each of the eight samples separately. These analyses exclude the two temporal samples and are based on 10,903 neutral SNPs. Values of r represent genetic autocorrelation coefficients, the dashed lines represent 95% CIs around the null hypothesis of randomly distributed genotypes, determined with 1000 permutations, and the error bars represent 95% CIs around r values, as determined by 1000 bootstraps. In the correlograms, the first x-intercept following significantly positive values of r (denoted by asterisks) indicates the extent of positive spatial autocorrelation or the genetic patch size

| Broad-scale structure
We uncovered two broad-scale patterns of population genetic structure in snapper across WA with our neutral SNP dataset.
Firstly, strong isolation by distance (IBD) was detected between the central west-coast of WA and the west-coast of South Australia (SA), indicating that dispersal is limited by coastal distance and probably primarily occurs in a stepping-stone fashion. Previous work with microsatellite markers and otolith microchemistry also uncovered IBD in snapper across WA (Edmonds et al., 1999;Fairclough et al., 2013;Gardner & Chaplin, 2011). IBD has also been detected in a number of other broadcast spawning marine species within our study range, including Roe's abalone (Hancock, 2000), WA dhufish (Berry et al., 2012) and saucer scallop (Kangas et al., 2019). However, that study lacked samples from locations between the Abrolhos Islands and Cockburn Sound so could not have detected a discontinuity just south of Lancelin. Additionally, although another otolith microchemistry study detected differentiation in snapper along the west coast of WA, a clear stock boundary was not detected between Cockburn Sound and Lancelin, potentially due to the often weak microchemistry differences in marine environments (Fairclough et al., 2013). Our study has therefore generated the clearest and most detailed account yet of the stock structure and connectivity of snapper in WA.
Generally, the stock discontinuities lie between locations within our study's range where the largest aggregations of spawning snapper have been observed -the Shark Bay region, Cockburn Sound (including the adjacent embayments Owen Anchorage and Warnbro Sound) and the SA gulfs (Spencer Gulf and Gulf St Vincent) -and therefore may indicate the importance of these spawning areas in maintaining spatially proximate snapper stocks. The discontinuities identified here delineate three broad-scale snapper stocks (i.e. the mid-west, south-west and south-coast stocks; see Figure 2). The size and spatial boundaries of the three snapper stocks are probably largely shaped by spawning timing and location, as well as oceanography and coastal topography. The mid-west stock is likely the larger of the two purely WA stocks because spawning grounds in the Shark Bay region are exposed to the Shark Bay Outflow (see Figure 1; Hetzel et al., 2018) and a strong Leeuwin Current (LC; Cresswell & Golding, 1980, Godfrey & Ridgway, 1985, facilitating the southward transport of pelagic life stages over potentially several hundreds of kilometres (Feng et al., 2010). In contrast, in the south-west, the magnitude of dispersal of PLS snapper from Cockburn Sound is probably limited because it is largely protected from shelf currents like the LC and Capes Current (see Figure 1; Pearce & Pattiaratchi, 1999). Instead, wind-driven gyres within the embayment (Steedman & Craig, 1983) facilitate the retention of eggs and larvae (Wakefield, 2010), limiting long-distance transport of PLS.
The boundary between the mid-west and south-west stocks lies between ecological bioregions (inshore: Southwest Shelf Transition and Southwest Shelf Province, offshore: Central Western Transition and Southwest Transition) delineated from patterns of species distributions and environmental features (Commonwealth of Australia, 2006). The boundary between these bioregions coincides with patterns in the distribution of invertebrates (Kott, 1952;O'Hara & Poore, 2000), fishes (Ayvazian & Hyndes, 1995;Hutchins, 1994;Last et al., 2011;Williams et al., 2001) and seaweeds (Wernberg et al., 2013). These studies, along with ours, indicate that there may be oceanographic features in the region capable of disrupting alongshore transport of PLS. Indeed, the marine environment around Lancelin is dynamic and comprises features that could achieve this.
South of the Abrolhos Islands ( Figure 1) the LC strengthens, leading to its instability and the production of eddies just north of Lancelin where the shelf narrows, causing the offshore movement of water (Feng et al., 2007;Pattiaratchi, 2006;Pearce & Griffiths, 1991).
Circulation patterns off Lancelin, along with the retention capability of the gyres in Cockburn Sound, may contribute to the stock boundary observed for snapper between the mid-west and south-west regions.
In contrast to the break between the mid-west and south-west stocks, a transition region occurred between the south-west and south-coast stocks. The heterogeneous Esperance sample occurred within this transition region, with approximately half of the sample having south-west ancestry and the other half having south-coast ancestry. The location of this transition region may in part reflect the magnitude of dispersal events from the SA gulfs (relative to those from further west). It is likely that this westward dispersal occurs post-recruitment, since thermal fronts form at the entrances of the SA gulfs during the spawning season, largely preventing gulf-shelf exchanges (Bruce & Short, 1990;Petrusevics, 1993;Vaz et al., 1990).
For example, in Nerita atramentosa (an intertidal snail with a pelagic larval duration of around four months), the percentage of larvae released during summer in the two SA gulfs that did not reach the boundary current ranged from 70% to 100% (Teske et al., 2015).
Hypothetically, long-distance dispersal from the SA gulfs may occur during the subadult stage, a conclusion of otolith microchemistry work (Fowler et al., 2005). We however acknowledge that the large sampling gap between Esperance and Ceduna (~1200 km) means that we cannot characterize the exact nature of the discontinuity between the south-west and south-coast stocks.

| Local-scale structure
The strength of an IBD relationship may be greater in some areas and weaker in others. The spatial autocorrelation analyses allowed us to investigate whether local-scale patterns of genetic structure in WA mirror the broad-scale signal of IBD. Our results indicated that the mid-west deviates from the broad-scale IBD pattern and constitutes a well-mixed stock of snapper, which is likely facilitated by the LC (flowing strongly during the winter, the mid-west spawning season). In contrast, in the south-west, we detected a pattern of spatial autocorrelation consistent with IBD, potentially facilitated by a weak LC during snapper's spring/summer spawning period in the region and the location of important spawning sites (i.e. protected embayments). Positive spatial autocorrelation occurred at distances up to ~300 km, suggesting that demographic connectivity between our sites on the western and southern coastlines of the region is probably limited and may occur in a stepping-stone fashion. King George Sound in Albany is an important spawning area for snapper along the south coast of WA (Wakefield et al., 2015). It is possible that this embayment, along with a number of smaller adjacent embayments, is most important for supplying parts of the southern coastline of the south-west (Neira & Potter, 1992;Potter et al., 1990). Overall, our spatial autocorrelation analyses suggest that the influence of oceanographic processes on the transport of PLS may be particularly important in shaping local-scale patterns of genetic structure in snapper in WA.

| Temporal variability in genetic structure
Analysis of the two temporally replicated samples uncovered differences in ancestry membership to the three identified stocks at Esperance, but not at Cockburn Sound. Temporal variation in population genetic structure has been reported for numerous marine species (Hogan et al., 2010;Jackson et al., 2018;Papetti et al., 2009;Quintero-Galvis et al., 2020;Watts et al., 2010).
When comparing the repeated Esperance samples, there was a significant reduction in ancestry characteristic of our SA sample over time and an increased contribution from stocks further west (i.e. mid-west and south-west). This shift suggests that the genetic structure of snapper off Esperance may largely depend on dynamics of snapper populations in SA and along the south coast of WA, including year-to-year recruitment strength. Indeed, interannual recruitment variability can be substantial in snapper and is considered to greatly influence population dynamics (Fowler & McGlennon, 2011;Francis, 1993;Hamer & Jenkins, 2004;Wakefield et al., 2015). Recent stock depletions and prolonged recruitment failure in the SA gulfs may account for the change in population genetic structure at Esperance (Fowler et al., 2020).
The shift in genetic composition of Esperance snapper also implies that the location of the genetic discontinuity in the south-east varies temporally. Overall, our temporal analyses indicate that incorporating samples collected at multiple timepoints in population genetic studies can reveal important biological information that may have otherwise gone undetected.

| Suitability of current stock boundaries based on population genomics
Employing a large genome-wide dataset gave us sufficient power to uncover patterns of population genetic structure and connectivity at multiple spatial scales. We can therefore confidently use our findings to inform the spatial assessment and management of snapper. and South Coast stocks (comprising Albany and Esperance; see Figure 1). For the West Coast, a whole of bioregion assessment is conducted along with two smaller area-based assessments, with the first including locations south of Lancelin and the second locations north of Lancelin to Kalbarri. The most recent assessment identified differences in stock status between these areas, related to variation in levels of exploitation (Fairclough et al., 2021). There are also known differences in biological characteristics (e.g. length and age) between the southern west coast and the central west coast (e.g. see Figure 1 inset). Consistent with this, we detected a stock boundary between Lancelin and Cockburn Sound, as well as a lack of genetic structure over an area extending ~800 km between the Gascoyne and Lancelin. Our results therefore suggest that managing the midwest stock separately may be beneficial. However, such an approach would need to consider the fact that the fisheries that operate in the region are multispecies and multisector and that other species may have differing population structuring.
In contrast to the mid-west, in the south-west of WA between Cockburn Sound and Albany, local-scale genetic structure was uncovered whereby genetic differentiation increased with coastal distance (i.e. an IBD pattern). Simulation work by Spies et al. (2015) showed that splitting an area characterized by IBD for management purposes leads to more favourable outcomes than if considered singularly (e.g. spawning biomass). Further, they demonstrated that the most optimal outcomes could be achieved by splitting a management region according to dispersal distance and fishing effort.
Others suggest that under a pattern of IBD, spatial scales of management might be based on the x-intercept of spatial autocorrelation correlograms, a value that indicates the distance at which individuals become genetically independent (Diniz-Filho & De Campos Telles, 2002;Östman et al., 2017). Considering our estimate of mean dispersal distance per generation determined from the IBD slope (~400 km), the x-intercept of our spatial autocorrelation correlogram (~300 km) and that fishing pressure is greater along the west coast near the Perth metropolitan region than the south coast (Gaughan & Santoro, 2021), the current management boundary on the south-west corner of the state (at 115.5°E, see Figure 1) appears suitable.
Along the south coast of WA, the region between the southwest corner and the SA border (i.e. the south-coast stock) is assessed and managed as a single unit (Newman et al., 2022). Here we detected significant gene flow across this region. However, the genetic composition of snapper in the eastern part of the region (i.e. at Esperance) varies temporally and therefore the location of the boundary between SA and WA snapper changes over time.
Considering this finding, and that catch data indicate that the biomass of snapper in the south east of WA is likely low relative to adjacent areas (Norriss et al., 2016), it may be appropriate to consider the eastern part of the southern coast of WA (i.e. east of Albany) separately for stock assessment and management purposes. The southern coast of WA presents a more complex situation for spatial assessment and management than the west coast and highlights the difficulties in translating complex biological processes into distinct units suitable for fisheries management. Population genomic work involving finer scale sampling along the southern coast of WA and the west coast of SA (i.e. between Albany and Ceduna) would be valuable for discerning the most suitable locations for drawing management boundaries. Additionally, our temporal analysis suggests that regular genetic monitoring would be beneficial to better understand the temporal and spatial dynamics of the stock boundaries uncovered here.

| CON CLUS IONS
The fluidity of the marine environment, along with the highly mobile and abundant nature of the species which occupy it, makes the characterization of evolutionary and demographic patterns in marine organisms difficult. Our study highlights the utility of large datasets of genetic markers in improving the understanding of spatial population structure and connectivity in marine species with high dispersal potential. Our findings likely relate to variation in physical factors, such as local ocean circulation and coastal geomorphology, as well as in biological factors, including timing of reproduction, spawning site preferences, recruitment variability and migratory behaviour.
Overall, our study demonstrates the value of population genomics in helping to improve the spatial management of fishery resources and therefore their long-term sustainability. Additionally, our work adds to the extensive body of literature showing that marine species are often characterized by complex population structure, rather than panmixia, and that population structure is not always static across time.

ACK N OWLED G M ENTS
This work was financially supported by the Australian Research Council (LP180100756 to LBB), Flinders University and by the PhD scholarships to AB provided by AJ and IM Naylon and the Playford Trust. We sincerely thank those who assisted with the collection of tissue samples, which included staff from the Department of Primary Industries and Regional Development WA and the South Australian Research and Development Institute. We also thank Michelle Gardner for providing the temporal sample from Esperance.
We thank Andrea Barceló and Diana-Elena Vornicu at the Molecular Ecology Lab Flinders University for their assistance in the laboratory.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.